Take-Home Exercise 1: Geospatial Analytics for Social Good: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar

1. Overview

1.1 Introduction

Since early 2021, Myanmar has been engulfed in a brutal civil war that erupted following the military coup on February 1, 2021, which ousted the democratically elected government led by Aung San Suu Kyi. The coup ignited widespread protests and a civil disobedience movement across the country, escalating into armed resistance as various ethnic armed groups and newly formed militias began confronting the military junta. The conflict has led to severe humanitarian crises, including widespread violence, displacement, and human rights abuses. The military’s response has been marked by extreme repression, including airstrikes and targeted attacks on civilian populations, exacerbating the suffering of ordinary Myanmar citizens. Despite international condemnation and calls for a return to democratic governance, the violence continues, deepening the political instability and socio-economic challenges faced by Myanmar.

In light of this ongoing crisis, I will be conducting a detailed geospatial analysis of the conflict events in Myanmar. This analysis aims to map and evaluate the spatial distribution and intensity of conflict incidents from 2021 to 2024. By utilizing geospatial techniques and spatial data, I will examine patterns of violence, identify hotspots of conflict, and assess the impact on civilian areas. This approach will provide valuable insights into the spatial dynamics of the conflict, helping to inform humanitarian responses and policy decisions.

1.2 Goal

We will be focusing on 4 event types: Battles, Explosion/Remote violence, Strategic developments, and Violence against civilians.

• Using appropriate function of sf and tidyverse packages, import and transform the downloaded armed conflict data and administrative boundary data into sf tibble data.frames.

• Using the geospatial data sets prepared, derive quarterly KDE layers.

• Using the geospatial data sets prepared, perform 2nd-Order Spatial Point Patterns Analysis.

• Using the geospatial data sets prepared, derive quarterly spatio-temporal KDE layers.

• Using the geospatial data sets prepared, perform 2nd-Order Spatio-temporal Point Patterns Analysis.

• Using appropriate tmap functions, display the KDE and Spatio-temporal KDE layers on openstreetmap of Myanmar.

• Describe the spatial patterns revealed by the KDE and Spatio-temporal KDE maps.1.1 The Data

1.3 Importing of Packages

Before we start off, we will have to import the necessary packages required for us to conduct our analysis.

We will be using the following packages:

• sf: Manages spatial vector data, enabling operations like spatial joins, buffering, and transformations for points, lines, and polygons.

• raster: Handles raster data, allowing for operations such as raster calculations, resampling, and visualization of spatial grids (e.g., elevation or satellite images).

• spNetwork: Analyzes spatial networks by modeling connectivity and movement within networks like road systems or utility grids.

• tmap: Creates and customizes thematic maps for spatial data visualization, including static and interactive maps with various map elements.

• tidyverse: A suite of packages for data manipulation (dplyr), visualization (ggplot2), and tidying (tidyr), facilitating a streamlined workflow for data analysis.

• RColorBrewer: Provides a range of color palettes for effective and aesthetically pleasing data visualization, including options for categorical, sequential, and diverging data.

• spatstat: Conducts spatial statistics and analysis, including point pattern analysis and spatial simulations, to study spatial distributions and interactions.

• sparr: Facilitates the analysis of spatial point patterns and spatial autoregressive models. It includes functions for fitting and analyzing spatial point processes, particularly useful for examining spatial dependencies and interactions in point pattern data.

• stpp : Provides tools for spatio-temporal point pattern analysis, including methods for modeling and analyzing the interactions and distributions of points over both space and time.

pacman::p_load(sf, raster, spNetwork, tmap, tidyverse, RColorBrewer, spatstat, sparr, stpp,lgcp)

1.4 The Data

For the purpose of this study, we will be using the following datasets. Particularly, I will be focusing on the quarterly armed conflict events from January 2021 until June 2024.

• Armed Conflict Location & Event Data of Myanmar between Jan 2021 to Jun 2024 from Armed Conflict Location & Event Data (ACLED), an independent, impartial, international non-profit organization collecting data on violent conflict and protest in all countries and territories in the world.

• Myanmar State and Region Boundaries with Sub-region from Myanmar Information Management Unit, MIMU

• OpenStreetMap dataset, which is an open-sourced geospatial dataset including shapefiles of important layers including road networks, forests, building footprints and many other points of interest.

• Myanmar District Boundaries from Myanmar Information Management Unit, MIMU

2. Data Wrangling

2.1 ACLED data

2.1.1 Importing Data

The below chunk of code is used to import ACLED conflict data from a CSV file and convert it into a geospatial data frame using longitude and latitude as coordinates.

In order to perform geoprocessing using two geospatial data, we need to ensure that both geospatial data are projected using similar coordinate system which is why in this case we will project it to WGS84 with the crs code of 32647 using st_transform.

Since the column event_date was stored as characters, dmy() is also used to format the event_date column into a standardized date format for further analysis.

st_as_sf is used to convert a data frame or other tabular data (like from CSVs, data frames, or tibbles) into a simple features (sf) object.

acled_sf <- read_csv("data/aspatial/2021-01-01-2024-06-30-Myanmar.csv") %>% 
  st_as_sf(coords = c(
    "longitude", "latitude"), crs = 4326) %>% 
  st_transform(crs= 32647)%>%
  mutate(event_date = dmy(event_date))

Now, let us check the CRS again by using the code chunk below.

st_crs(acled_sf)

Coordinate Reference System:
  User input: EPSG:32647 
  wkt:
PROJCRS["WGS 84 / UTM zone 47N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 47N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",99,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
        BBOX[0,96,84,102]],
    ID["EPSG",32647]]

2.1.2 Data Pruning and Structuring

2.1.2.1 Extracting Data for the 4 Main Event Types

Since the data we have encapsulates all 6 types of events and we are only focusing on the 4 event types that I have indicated above, we will only be extracting the data related to that.

acled_sf <- acled_sf %>%
  filter(event_type %in% c("Battles", "Strategic developments", "Violence against civilians", "Explosions/Remote violence"))

Using unique(), we can check that only rows that fall under the 4 event types have been kept.

unique(acled_sf$event_type)

[1] "Battles"                    "Strategic developments"    
[3] "Violence against civilians" "Explosions/Remote violence"

2.1.2.2 Data Transformation for Quarterly Analysis

Since we will be analyzing it by quarters, let’s use mutate() to extract the information. Using quarter() from the lubridate package, we can derive the quarters from the event_date.

acled_sf <- acled_sf %>% mutate(quarter = quarter(event_date))

Using colnames() we can see that the new column quarter is added.

colnames(acled_sf) 

 [1] "event_id_cnty"      "event_date"         "year"              
 [4] "time_precision"     "disorder_type"      "event_type"        
 [7] "sub_event_type"     "actor1"             "assoc_actor_1"     
[10] "inter1"             "actor2"             "assoc_actor_2"     
[13] "inter2"             "interaction"        "civilian_targeting"
[16] "iso"                "region"             "country"           
[19] "admin1"             "admin2"             "admin3"            
[22] "location"           "geo_precision"      "source"            
[25] "source_scale"       "notes"              "fatalities"        
[28] "tags"               "timestamp"          "geometry"          
[31] "quarter"           

2.1.2.3 Removal of Redundant Columns

To enhance the efficiency of our dataset, we will also remove redundant columns. This practice reduces memory usage and processing time while simplifying the analysis. By focusing on only the relevant data, we minimize complexity and ensure clearer, more focused results. Examples of the columns we will be removing are: time_precision, inter1 etc. With this, we are able to reduce to a total of 20 columns.

# Define columns to be removed
columns_to_remove <- c("time_precision", "inter1", "inter2", "iso", "source", "source_scale", "notes", "tags", "region", "geo_precision", "source_scale","country")

# Remove columns only if they exist in the dataframe
acled_sf <- acled_sf %>%
  dplyr::select(-all_of(columns_to_remove[columns_to_remove %in% names(acled_sf)]))

This is an overview of all the columns left:

colnames(acled_sf)

 [1] "event_id_cnty"      "event_date"         "year"              
 [4] "disorder_type"      "event_type"         "sub_event_type"    
 [7] "actor1"             "assoc_actor_1"      "actor2"            
[10] "assoc_actor_2"      "interaction"        "civilian_targeting"
[13] "admin1"             "admin2"             "admin3"            
[16] "location"           "fatalities"         "timestamp"         
[19] "geometry"           "quarter"           

2.1.2.4 Grouping of Data By Year, Quarter, Event Type

This chunk of code groups the ACLED Data by quarter, event and year, allowing us to better manage and access the data later on.

grped_acled_sf <- acled_sf %>%
  group_by(year, quarter, event_type) %>%
  summarize(event_count = n(), .groups = 'drop')

2.2 Myanmar Boundary and Sub-Region Dataset from MIMU

The code chunk below uses st_read() function of sf package to import mmr_polbnda2_adm1_250k_mimu_1 shapefile into R as a polygon feature data frame.

mmrsr_sf <- st_read(dsn="data/geospatial/mmr_polbnda2_adm1_250k_mimu_1/", 
                   layer="mmr_polbnda2_adm1_250k_mimu_1")

Reading layer `mmr_polbnda2_adm1_250k_mimu_1' from data source 
  `/Users/georgiaxng/georgiaxng/is415-handson/Take-home_Ex/Take-home_Ex01/data/geospatial/mmr_polbnda2_adm1_250k_mimu_1' 
  using driver `ESRI Shapefile'
Simple feature collection with 18 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84

As shown below, currently the data is in geographic coordinate system (latitude/longitude), we will need to transform it to a projected CRS before proceeding.

st_crs(mmrsr_sf)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["geodetic longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]

mmrsr_sf <- st_transform(mmrsr_sf, crs = 32647)
st_crs(mmrsr_sf)

Coordinate Reference System:
  User input: EPSG:32647 
  wkt:
PROJCRS["WGS 84 / UTM zone 47N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 47N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",99,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 96°E and 102°E, northern hemisphere between equator and 84°N, onshore and offshore. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Myanmar (Burma). Russian Federation. Thailand."],
        BBOX[0,96,84,102]],
    ID["EPSG",32647]]

To check that all the spatial objects are valid, st_is_valid() is utilized. As shown below, all of the objects are valid.

st_is_valid(mmrsr_sf)

 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[16] TRUE TRUE TRUE

This is a brief overview of how the map of Myanmar looks like.

plot(mmrsr_sf)

2.3 Mapping the Geospatial data sets

Now, let’s plot the spatial map to gain an initial understanding of the geographic distribution of all of the 4 conflict events by year. As observed, with such a large volume of data, identifying patterns through visual inspection alone can be challenging. This highlights the importance of conducting further analysis to uncover trends and achieve a deeper understanding of the spatial dynamics of these events.

From a brief examination of the map, it’s evident that the region south of Sagaing exhibits the highest concentration of conflicts overall. This area shows a notably higher density of conflict events compared to other regions, marking it as a significant hotspot. However, to fully grasp these patterns and their implications, additional analysis is needed.

2021

# Filter data for 2021
acled_2021 <- acled_sf %>% filter(year == 2021)

# Plot map for 2021
tm_shape(mmrsr_sf) +
  tm_polygons() + 
  tm_shape(acled_2021) + 
  tm_dots() +
  tm_layout(main.title = "ACLED Conflict Events in 2021", main.title.size = 1.2)

2022

# Filter data for 2022
acled_2022 <- acled_sf %>% filter(year == 2022)

# Plot map for 2022
tm_shape(mmrsr_sf) +
  tm_polygons() + 
  tm_shape(acled_2022) + 
  tm_dots() +
  tm_layout(main.title = "ACLED Conflict Events in 2022", main.title.size = 1.2)

2023

# Filter data for 2023
acled_2023 <- acled_sf %>% filter(year == 2023)

# Plot map for 2023
tm_shape(mmrsr_sf) +
  tm_polygons() + 
  tm_shape(acled_2023) + 
  tm_dots() +
  tm_layout(main.title = "ACLED Conflict Events in 2023", main.title.size = 1.2)

2024

# Filter data for 2024
acled_2024 <- acled_sf %>% filter(year == 2024)

# Plot map for 2024
tm_shape(mmrsr_sf) +
  tm_polygons() + 
  tm_shape(acled_2024) + 
  tm_dots() +
  tm_layout(main.title = "ACLED Conflict Events in 2024", main.title.size = 1.2)

2.4 Myanmar Sub-regions: For Reference and Context

For reference to the different sub-regions in Myanmar, we can refer to the labelled map plotted below:

num_colors <- length(unique(mmrsr_sf$ST))
colors <- brewer.pal(n = num_colors, name = "Set3")

tm_shape(mmrsr_sf) +
  tm_polygons(col = "ST", palette = colors) +  # Apply color palette to polygons
  tm_text("ST", size = 1, col = "black", bg.color = "white", just = c("center", "center"),  xmod = 0, ymod = 0) + tm_layout(main.title = "Myanmar",
            main.title.position = "center",
            main.title.size = 1.6,
            legend.outside = TRUE,
            frame = TRUE)+
    tm_legend(title = "Sub-regions")  # Set custom legend title

2.5 Myanmar District Data

This chunk of code imports the geospatial data for the districts of Myanmar from MIMU.

mmrsr_district_sf <- st_read(dsn="data/geospatial/mmr_district_mimu/", 
                   layer="mmr_polbnda_adm2_250k_mimu")

Reading layer `mmr_polbnda_adm2_250k_mimu' from data source 
  `/Users/georgiaxng/georgiaxng/is415-handson/Take-home_Ex/Take-home_Ex01/data/geospatial/mmr_district_mimu' 
  using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84

mmrsr_district_sf <- st_transform(mmrsr_district_sf, crs = 32647)

plot(mmrsr_district_sf)

2.6 Setting Seed

To ensure the reproducibility, we will be using a fixed seed for our randomisation.

set.seed(123456)

3. Kernel Density Estimation

To analyze and find the trends in the spatial distribution of conflict events in Myanmar, we will be using Kernel Density Estimation (KDE). KDE is a non-parametric way to estimate the probability density function of a random variable, making it particularly useful in spatial analysis for identifying areas with high concentrations of events.

To determine the optimal bandwidth and kernel parameters for KDE, we will be using a sample of the data. By focusing on a smaller subset, such as specific event types or time periods, we can efficiently explore different parameter combinations and choose the one that provides the most meaningful results before scaling the analysis to the entire dataset.

3.1 Creating owin object

To improve the accuracy of our spatial analysis, we will exclude the smaller islands that are too minor to contribute meaningfully to the Kernel Density Estimation (KDE). We will focus on the main island by merging the geometry and selecting the relevant polygon.

In this code snippet, we combine all intersecting geometries in the mmrsr_sf object into a single polygon using st_union(). We then use st_cast() to ensure that the result is a polygon type. From this merged polygon, we select the row with the highest number of features, which typically corresponds to the main island or the most prominent feature in the dataset.

This selected polygon is then converted to a window object using as.owin(), which allows us to use it in spatial analyses, such as Kernel Density Estimation (KDE).

merged_mmr <- st_union(mmrsr_sf) %>%
    st_cast("POLYGON")
merged_mmr <- merged_mmr[c(9)]

mmr_owin <- as.owin(merged_mmr)

Here is the plot of the main island, without its outer islands.

plot(mmr_owin)

3.2 Selection of Sample Data

To choose a representative sample, we will first calculate the median number of data points across different groupings, such as year, quarter, and event type. This will help us understand the typical size of our data subsets. Using this median value, we will identify the subset of data with an event_count closest to this median.

By selecting a sample that aligns with this typical value, we ensure that our sample accurately reflects the general distribution of the dataset. We will then use this representative sample to determine the optimal parameters for Kernel Density Estimation (KDE), such as bandwidth and kernel type. This approach will help us achieve more reliable and meaningful results in our spatial analysis.

Here is the R code to select the representative sample. In this case, the resulting sample data we have arrived at is for “Battles” in the second quarter of 2023.

# Calculate mean and median of event counts
median_event_count <- median(grped_acled_sf$event_count)

# Find the row with event_count closest to the median
sample_acled_sf <- grped_acled_sf %>%
  mutate(distance_to_median = abs(event_count - median_event_count)) %>%
  arrange(distance_to_median) %>%
  slice(1)

# Print the closest row(s)
sample_acled_sf

Simple feature collection with 1 feature and 5 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -174035.7 ymin: 1103500 xmax: 518300.4 ymax: 3006373
Projected CRS: WGS 84 / UTM zone 47N
# A tibble: 1 × 6
   year quarter event_type event_count                                  geometry
  <dbl>   <int> <chr>            <int>                          <MULTIPOINT [m]>
1  2023       2 Battles            814 ((-174035.7 2284958), (-152179.1 2277258…
# ℹ 1 more variable: distance_to_median <dbl>

Here, we converted the representative sample data to ppp format and checked for duplicates. In this case, no duplicates were found, confirming the integrity of our sample for further KDE analysis.

sample_acled_ppp <- as.ppp(st_coordinates(sample_acled_sf),st_bbox(sample_acled_sf))
any(duplicated(sample_acled_ppp))

[1] FALSE

sample_acled_ppp <- sample_acled_ppp[mmr_owin]

An overview of the point pattern plot for the sample data.

plot(sample_acled_ppp, pch = 16, cex = 0.8) 

3.3 Determining the Optimal Bandwidth & Sigma

3.3.1 Understanding Sigma & Kernel in Kernel Density Estimation (KDE)

Kernel:

The kernel in Kernel Density Estimation (KDE) is a smoothing function that estimates the probability density function of a dataset. The choice of kernel function affects how the smoothing is applied to each data point. Common kernels include:

• Gaussian Kernel: Assumes a normal distribution around each data point, providing a smooth, continuous density estimate.

• Epanechnikov Kernel: Parabolic in shape, it minimizes mean squared error but is less smooth than Gaussian.

• Uniform Kernel: Assumes constant density within a fixed distance around each data point, simpler but less smooth.

• Triangular Kernel: Offers a clear density representation with sharp edges but less smoothness.

The choice of kernel affects the density estimate’s smoothness and shape, with the Gaussian kernel being a versatile and commonly used option.

Sigma:

The sigma argument in the density() function represents the bandwidth of the kernel in KDE. It controls the level of smoothing applied to the density estimate:

• What Sigma Represents: Sigma defines the standard deviation of the kernel function, determining the width of the “smoothing window” around each data point.

• Effect on KDE:

  • Smaller Sigma: Leads to less smoothing, producing a more sensitive and detailed density plot, but can highlight noise.

  • Larger Sigma: Results in greater smoothing, providing a broader view of density but might obscure finer details.

• Choosing Sigma: An appropriate sigma balances detail and smoothness. Methods like bw.diggle, bw.ppl, and bw.scott can help determine the optimal sigma value.

In summary, the kernel determines the density estimate’s shape, while sigma controls the smoothing level. Both are crucial for accurately reflecting the spatial distribution of data, and selecting the right parameters is essential for meaningful KDE results.

3.3.2 Computing Default Kernel Density Estimation

First, since the original data is in meters, we need to rescale it to kilometers to facilitate more meaningful spatial analysis and visualization. By converting the units to kilometers, we ensure that the scale of the Kernel Density Estimation (KDE) and subsequent analyses align with the geographical extent of the study area.

sample_acled_ppp.km <- rescale.ppp(sample_acled_ppp, 1000, "km")

To visualize the spatial distribution of conflict events, here is the Kernel Density Estimation (KDE) plot generated using the default settings, as shown in the code. While default settings can be sufficient for some analyses, the resulting KDE plot here appears oversmoothed, which may obscure finer details and potential small-scale clusters.

This oversmoothing effect reduces the ability to detect subtle patterns in the conflict events, highlighting the importance of adjusting parameters like bandwidth (sigma) to better capture the underlying spatial structure. Thus, to improve the accuracy and granularity of the KDE, we will need to continue experimenting with different bandwidth values and kernel functions.

par(mar = c(0,1,1,1))
kde_default_destination <- density(sample_acled_ppp.km)
plot(kde_default_destination,main = "Default Density KDE for ACLED Myanmar 2023 Q2 'Battles'")

3.4 Experimenting With Fixed Bandwidth

To effectively choose the appropriate bandwidth for Kernel Density Estimation (KDE), several automatic bandwidth selection methods can be employed. Each method has its own approach and focus, which can significantly impact the resulting density estimate. Below are brief descriptions of four common bandwidth selection methods:

• bw.CvL: Calculates bandwidth using cross-validation to minimize the integrated mean squared error (IMSE) of the density estimate. This method aims to balance detail and smoothness by optimizing the bandwidth based on a global error metric, which might lead to a broad smoothing effect.

• bw.scott: Applies a rule-of-thumb bandwidth based on the sample size and dimensionality of the data. This method is often effective for larger datasets, providing a more generalized estimate by assuming isotropic smoothing (uniform smoothing in all directions).

• bw.ppl: Uses a bandwidth tailored for point pattern analysis, considering the spatial distribution of data points. This approach focuses on capturing local variations and spatial patterns, resulting in a more localized smoothing effect compared to other methods.

• bw.diggle: Employs a method specifically designed for point patterns to minimize the variance of the density estimate. This results in a smaller bandwidth, leading to a finer and more detailed density plot that may capture small-scale variations but could also emphasize noise.

3.4.1 Automatic Bandwidth Methods

Lets begin by examining the bandwidth values obtained from various automatic bandwidth selection methods using the following code snippet.

bw.ppl(sample_acled_ppp.km)

   sigma 
21.20514 

bw.diggle(sample_acled_ppp.km)

   sigma 
8.459462 

bw.CvL(sample_acled_ppp.km)

   sigma 
90.92206 

bw.scott(sample_acled_ppp.km)

  sigma.x   sigma.y 
 51.29856 133.58599 

Based on the sigma values obtained from the automatic bandwidth calculation methods, we can draw the following inferences about the Kernel Density Estimation (KDE) for the sample data:

• bw.ppl: sigma = 21.21
  A smaller sigma value suggests a more localized smoothing effect. This method is designed for point pattern analysis and aims to capture spatial patterns with greater precision, highlighting finer details within the data.

• bw.diggle: sigma = 8.46
  As the smallest sigma value, it implies a very narrow smoothing window. This method focuses on minimizing variance and tends to produce a detailed density plot. However, it may also accentuate noise or small fluctuations in the data.

• bw.CvL: sigma = 90.92
  This large value indicates a broad smoothing window. Consequently, each data point affects a wide surrounding area, leading to a smoother KDE that may obscure fine details and small-scale variations.

• bw.scott: sigma.x = 51.30, sigma.y = 133.59
  The method provides different bandwidths for the x and y dimensions, suggesting anisotropic smoothing. The significant difference between these values indicates varying spatial variations in different directions, with a broader smoothing effect along the y-direction.

3.4.2 Plotting KDE for Different Bandwidth & Kernel

To better understand these differences, we will plot the KDE results using each bandwidth value and different kernels. This visual comparison will help us assess how each bandwidth setting and kernel type influences the density estimates, allowing us to determine which combination provides the most meaningful representation of our data.

In practice, choosing the optimal KDE bandwidth is not straightforward, as there is no one-size-fits-all approach. Many studies emphasize that determining the best bandwidth often relies on visually comparing various settings and kernel types to assess their effectiveness and select the most appropriate one for the specific data at hand.

3.4.2.1 Overview of KDE Maps with Gaussian Kernel and Various Bandwidths

Here is an overview of how in general the KDE maps generated with the respective bandwidths looks like using the Gaussian kernel.

kde_diggle <- density(sample_acled_ppp.km, bw.diggle(sample_acled_ppp.km))
kde_CvL <- density(sample_acled_ppp.km, bw.CvL(sample_acled_ppp.km))
kde_scott <- density(sample_acled_ppp.km, bw.scott(sample_acled_ppp.km))
kde_ppl <- density(sample_acled_ppp.km, bw.ppl(sample_acled_ppp.km))

par(mar = c(1,1,1,1.5),mfrow = c(2,2))
plot(kde_diggle,main = "kde_diggle")
plot(kde_CvL,main = "kde_CvL")
plot(kde_scott,main = "kde_scott")
plot(kde_ppl,main = "kde_ppl")

3.4.2.2 Comparing KDE Results Across Different Bandwidths and Kernels

The code chunk below plots different KDE plots for a various different pairings of kernels and bandwidths.

bw.ppl

par(mfrow=c(2,2))
plot(density(sample_acled_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="gaussian"), main="Gaussian")
plot(density(sample_acled_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="epanechnikov"), main="Epanechnikov")
plot(density(sample_acled_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="quartic"), main="Quartic")
plot(density(sample_acled_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="disc"), main="Disc")

bw.diggle

par(mfrow=c(2,2))
plot(density(sample_acled_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian"), main="Gaussian")
plot(density(sample_acled_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="epanechnikov"), main="Epanechnikov")
plot(density(sample_acled_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="quartic"), main="Quartic")
plot(density(sample_acled_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="disc"), main="Disc")

bw.CvL

par(mfrow=c(2,2))
plot(density(sample_acled_ppp.km, sigma=bw.CvL, edge=TRUE, kernel="gaussian"), main="Gaussian")
plot(density(sample_acled_ppp.km, sigma=bw.CvL, edge=TRUE, kernel="epanechnikov"), main="Epanechnikov")
plot(density(sample_acled_ppp.km, sigma=bw.CvL, edge=TRUE, kernel="quartic"), main="Quartic")
plot(density(sample_acled_ppp.km, sigma=bw.CvL, edge=TRUE, kernel="disc"), main="Disc")

bw.scott

par(mfrow=c(2,2))
plot(density(sample_acled_ppp.km, sigma=bw.scott, edge=TRUE, kernel="gaussian"), main="Gaussian")
plot(density(sample_acled_ppp.km, sigma=bw.scott, edge=TRUE, kernel="epanechnikov"), main="Epanechnikov")
plot(density(sample_acled_ppp.km, sigma=bw.scott, edge=TRUE, kernel="quartic"), main="Quartic")
plot(density(sample_acled_ppp.km, sigma=bw.scott, edge=TRUE, kernel="disc"), main="Disc")

After evaluating the KDE plots, I opted for the Gaussian kernel with the bandwidth provided by bw.ppl because this method provided a more localized smoothing effect, better capturing the spatial patterns in the data. In comparison, the bw.CvL method yielded a much larger bandwidth, leading to excessive smoothing, while bw.scott showed considerable disparity between x and y dimensions, suggesting anisotropic smoothing. The bw.diggle method, though precise, resulted in very narrow smoothing that might overemphasize noise.

Despite this, I found that the Gaussian KDE with bw.ppl still appeared to be undersmoothed. Consequently, further optimization is necessary to refine the bandwidth setting and achieve a more balanced density estimate that accurately represents the underlying spatial patterns in the data.

---

3.4.3 Fine-Tuning KDE Parameters

Since the Gaussian KDE using bw.ppl seems to be undersmoothed, we will need to adjust the sigma to achieve the desired results. To address this undersmoothing issue, we can increase the sigma value to widen the smoothing window, as demonstrated below. Through continued experimentation, we have determined these final parameters.

adjusted_bw <- bw.ppl(sample_acled_ppp.km) * 1.5
kde_adjusted <- density(sample_acled_ppp.km, sigma=adjusted_bw, edge=TRUE, kernel="gaussian")
plot(kde_adjusted)

Here is a comparison between the point pattern and the KDE plot of the sample data, this comparison allows us to evaluate how effectively the KDE captures the spatial distribution of the events and highlights areas where the KDE might need further refinement.

par(mar = c(1,0,1,0))
par(mfrow=c(1,2))
plot(sample_acled_ppp, pch = 16, cex = 0.5) 
plot(kde_adjusted)

Value of the adjusted bandwidth:

adjusted_bw

   sigma 
31.80771 

4. Quarterly KDE Layers

4.1 2021

4.1.1 Q1

This chunk of code generates the KDE layers using the adjusted bandwidth we have previously counted and the gaussian kernel, plotting the resulting KDE layers by each quarter.

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2021
  current_quarter = 1
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

Note

• Battles: Distinct hotspot in Shan (North), which indicates this is likely the main battlefield during this quarter.

• Explosions/Remote Violence: Appears to have a wider spatial distribution, with more intense activity seen along the southern, western and central-eastern parts of the country. This suggests that attacks carried out here were more widespread and not confined to one region. This could imply different actors using remote violence across various territories or a strategic shift towards this type of warfare in multiple conflict zones.

• Strategic Developments, Violence Against Civilians: We can observe that both events display a similar KDE map, indicating that both events could have been carried out concurrently, hence leading to the similarity. The co-occurrence of strategic developments and civilian violence is concerning and could imply that the areas experiencing strategic military activity are also sites where civilians are more vulnerable to harm. These could be areas under military occupation, zones of forced displacement, or regions where there is a breakdown of law and order.

4.1.2 Q2

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2021
  current_quarter = 2
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

Note

* Battles: Presents a more spatially distributed distribution, with new hotspots at the Chin/Sagaing/Magway intersection.

* Explosions/Remote Violence: The intensity of explosions and remote violence has shifted southward, with the most significant hotspot now in Yangon, the capital of the country. The shift to Yangon could also be a sign of more strategic attacks aimed at disrupting economic or political centers of power, following the military coup in early 2021.

* Strategic Developments, Violence Against Civilians: Both events still present similar spatial distribution, reinforcing the idea that military strategic moves are closely linked to targeted violence against civilians.

4.1.3 Q3

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2021
  current_quarter = 3
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

Note

All 4 events shows a distinct hotspot at southern Sagaing, which appears to be the main battlefield and where most of the conflicts happen for 2021 Q4. Similarly, Yangon appears to have a significant amount of conflicts happening as well.

4.1.4 Q4

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2021
  current_quarter = 4
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.1.5 Overall Insights for 2021

Fluid Conflict Zones: The hotspots of conflict in Myanmar were more fluid and dynamic, with notable shifts in violence across various regions throughout the year. This shifting pattern likely reflects the early stages of the Myanmar civil conflict, which began in full force following the February 2021 military coup. There is a strong focus placed on Yangon and Southern Sagaing in particular.

4.2 2022

4.2.1 Q1

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2022
  current_quarter = 1
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.2.2 Q2

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2022
  current_quarter = 2
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.2.3 Q3

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2022
  current_quarter = 3
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.2.4 Q4

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2022
  current_quarter = 4
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.2.5 Overall Insights for 2022

Stabilization of Conflict Zones: The KDE plot for all quarters for the year of 2022 appears to stay identical, with majority of the conflicts being in the vicinity of southern Sagaing. This consistency indicates a stabilization of the conflict, with fewer fluctuations in the geography of violence. This suggests that the main battleground has become more fixed, with both government forces and resistance groups focusing their efforts in this region.

Implications: The entrenchment of conflict in southern Sagaing also likely reflects a combination of military control over key areas and entrenched resistance by local groups. This stabilization might also indicate a stalemate or protracted conflict, where neither side is able to decisively push the front lines or gain control of new territories.

4.3 2023

4.3.1 Q1

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2023
  current_quarter = 1
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.3.2 Q2

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2023
  current_quarter = 2
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.3.3 Q3

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2023
  current_quarter = 3
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.3.4 Q4

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2023
  current_quarter = 4
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

4.3.5 Overall Insights for 2023

Consistent hotspots: There’s a persistent area of high activity in the region around Southern Sagaing across all event types and quarters, similar to the previous year.

Event type variations: Different event types show slightly different patterns. For example, “Strategic developments” tend to be more concentrated, while “Battles” and “Explosions/Remote violence” show more widespread activity.

Geographical spread: Activity is not uniform across the country. The southern peninsular region generally shows less activity, while the central and northern regions show more.

Q4 escalation: There appears to be an increase in activity across all event types in Q4, particularly in the northern regions.

4.4 2024

4.4.1 Q1

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2024
  current_quarter = 1
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

Note

Less conflicts: Conflict events appears to have lessen significantly, as evidenced from the decrease in the intensity level of the hotspots across Myanmar. This decline is particularly apparent for “Explosions/Remote Violence,” where the once prominent hotspots have become less concentrated and more dispersed compared to previous quarters.

4.4.2 Q2

  par(mar = c(1,0,1,0))
  par(mfrow=c(2,2))
  current_year = 2024
  current_quarter = 2
  quarter_data <- grped_acled_sf %>% filter(year == current_year, quarter == current_quarter) 
  # Loop through each event
  for (i in seq_len(nrow(quarter_data))) {
    
    filtered_sf <- quarter_data[i, ]

    # Extract the current year, quarter, and event type
    current_event <- filtered_sf$event_type
    
    # Convert to ppp object
    filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
    
    if(any(duplicated(filtered_ppp))){
      rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
    }
    filtered_ppp <- filtered_ppp[mmr_owin]
    # Rescale to kilometres
    filtered_ppp.km <- rescale.ppp(filtered_ppp, 1000, "km")
  
    # Conduct KDE with the chosen bandwidth
    kde_result <- density(filtered_ppp.km, sigma = adjusted_bw, edge = TRUE, kernel = "gaussian")
    
    # Plot the KDE
    plot(kde_result, main = paste(current_event," ",current_year, " Q", current_quarter, sep = ""))
}

Note

Battles: Hotspot in southern Sagaing has decreased in intensity

Violence Against Civilians: Hotspot in Sagaing remains, new hotspots are also appearing on the southern part of Myanmar.

5. 2nd-Order Spatial Point Patterns Analysis

In this section, we will conduct second-order spatial point pattern analysis (PPA), which examines interactions between events at various spatial scales. By using statistical functions such as G, F, K, and L, we can quantify and explore clustering patterns, distances between events, and the spatial distribution of points relative to a given area. These functions help reveal whether the events exhibit random, clustered, or regular patterns, providing deeper insights into the underlying spatial structure of conflict events in Myanmar.

Here’s a detailed description of the second-order spatial point pattern analysis functions:

• G-function (Nearest Neighbor Distribution Function): The G-function describes the cumulative distribution of the distances from each event to its nearest neighboring event. It highlights how clustered events are by comparing observed distances with those expected from a random spatial pattern.

• F-function (Empty Space Function): The F-function examines the distribution of distances from arbitrary (random) locations to the nearest event. It helps in evaluating how well events fill the study area. A large gap between random points and events indicates spatial heterogeneity or clustering.

• K-function (Ripley’s K-function): This function estimates the number of events within a certain distance of each other, relative to what would be expected under complete spatial randomness (CSR). It accumulates the number of events within varying distances, thus detecting clustering or dispersion at multiple spatial scales. If the observed K-function is higher than expected under CSR, it indicates clustering.

• L-function (Besag’s L-function): The L-function is a variance-stabilized transformation of the K-function. It linearizes the expected K-values, making it easier to interpret deviations from randomness. An L-function above the CSR line indicates clustering, while values below suggest dispersion.

5.1 Analysis Background

For the analysis for the whole of Myanmar, I will be primarily using the G-function to conduct second-order spatial point pattern analysis. I chose it for its computational efficiency, straightforward interpretation of nearest-neighbor distances, and suitability for analyzing a single type of event. The K and L-functions, while valuable, are more resource-intensive, making them not ideal for my large dataset, while the F-function is less applicable to our single-type event analysis.

For the analysis of Sagaing however, I will be using the K-function to conduct the relevant analysis.

5.1.1 Monte Carlo Simulation

I will also be conducting a Monte Carlo simulation for testing Complete Spatial Randomness (CSR) in the spatial distribution of conflict events. This involves generating multiple random point patterns under CSR to create a range of expected spatial statistics. By comparing the observed G-function of the conflict events against these simulated envelopes, I can determine if the observed clustering or dispersion is significantly different from what would be expected by chance.

The envelope function provides upper and lower bounds for the G-function, helping to assess whether the spatial pattern of conflict events is significantly clustered or dispersed compared to random patterns. This method ensures a rigorous evaluation of spatial randomness and enhances the interpretation of observed patterns.

5.1.2 Example of G-Function Plot

This code snippet below provides an example of a G-function plot of a sample of the data, calculated using an envelope derived from 99 simulations. We will also be keeping the confidence level of the envelopes to the default 95%.

This approach is used to evaluate the spatial distribution of events and assess how the observed point pattern compares to what might be expected under complete spatial randomness. Note that we will be using a set seed to ensure the reproducibility of our code.

set.seed(123456)
test_second_sf <- grped_acled_sf %>% filter(event_type == "Battles", year == 2023, quarter ==2)
test_second_ppp<-as.ppp(st_coordinates(test_second_sf), st_bbox(test_second_sf))

if(any(duplicated(test_second_ppp))){
  test_second_ppp <- rjitter(test_second_ppp, retry=TRUE, nsim=1, drop=TRUE)
}

test_second_ppp <- test_second_ppp[mmr_owin]
test_second.csr <- envelope(test_second_ppp, Gest, nsim = 99)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
99.

Done.

plot(test_second.csr, xlim=c(0,10000), main="Example of G-Function Plot")

Tip

The plot typically shows several key elements related to both the observed data and the simulated patterns:

• G^obs (Observed G-function):This represents the empirical G-function calculated from your actual point pattern data. It shows the cumulative distribution of the nearest-neighbor distances between points in the dataset.

• G^theo (Theoretical G-function): Theoretical G-function that would be expected if the point pattern follows Complete Spatial Randomness (CSR). Under CSR, the points are distributed independently and uniformly across the study area.

• G^hi (Upper Envelope):

  • Upper bound of the G-function values generated from the random simulations (based on nsim, the number of simulations). It defines the maximum value that the G-function could take under CSR in these simulations.

  • If the observed G-function (G^obs) exceeds G^hi, it suggests clustering—points are closer together than expected by chance.

• G^lo (Lower Envelope):

  • Lower bound of the G-function values from the random simulations. It defines the minimum G-function value expected under CSR.

  • If the observed G-function (G^obs) falls below G^lo, it suggests dispersion or regularity—points are more evenly spaced than expected by chance.

• Shaded area: Represents the 95% confidence interval

5.2 Whole of Myanmar

5.2.1 Hypothesis

To confirm the observed spatial patterns, I will conduct a hypothesis test using the G-function along with Monte Carlo simulation envelopes. The test will be designed as follows:

• H₀: The distribution of conflict events in Myanmar is randomly distributed (follows Complete Spatial Randomness, CSR).

• H₁: The distribution of conflict events in Myanmar is not randomly distributed.

• If the observed G-function lies outside the envelope, we will reject the null hypothesis, suggesting the presence of spatial clustering or dispersion.

• A 95% confidence interval will be used for the envelopes.

5.2.2 Battles

This chunk of code generates a G-function plot for all quarters of a particular event, using an envelope derived from 99 simulations for significance testing. Then, it saves the results to a list to plot it later on.

set.seed(123456)
battle_data_sf <- grped_acled_sf %>% filter(event_type == "Battles")

# Initialize an empty list to store the results
battle_G_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(battle_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- battle_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[mmr_owin]
  G_CK.csr <- envelope(filtered_ppp, Gest, nsim = 99)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  battle_G_CK_list[[name]] <- G_CK.csr
}

write_rds(battle_G_CK_list, "data/rds/battle_G_CK_list.rds")

This chunk of code plots the G-function plot for all quarters.

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(battle_G_CK_list)) {
  g_ck <- battle_G_CK_list[[i]]
  plot(g_ck, xlim=c(0,10000), main=names(battle_G_CK_list)[i])
}

5.2.3 Explosions/Remote Violence

set.seed(123456)
explosions_data_sf <- grped_acled_sf %>% filter(event_type == "Explosions/Remote violence")

# Initialize an empty list to store the results
explosions_G_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(explosions_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- explosions_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[mmr_owin]
  G_CK.csr <- envelope(filtered_ppp, Gest, nsim = 99)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  explosions_G_CK_list[[name]] <- G_CK.csr
}
write_rds(explosions_G_CK_list, "data/rds/explosions_G_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(explosions_G_CK_list)) {
  g_ck <- explosions_G_CK_list[[i]]
  plot(g_ck, xlim=c(0,10000), main=names(explosions_G_CK_list)[i])
}

5.2.4 Strategic Developments

set.seed(123456)
strategic_data_sf <- grped_acled_sf %>% filter(event_type == "Strategic developments")

# Initialize an empty list to store the results
strategic_G_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(strategic_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- strategic_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[mmr_owin]
  G_CK.csr <- envelope(filtered_ppp, Gest, nsim = 99)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  strategic_G_CK_list[[name]] <- G_CK.csr
}
write_rds(strategic_G_CK_list, "data/rds/strategic_G_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(strategic_G_CK_list)) {
  g_ck <- strategic_G_CK_list[[i]]
  plot(g_ck, xlim=c(0,10000), main=names(strategic_G_CK_list)[i])
}

5.2.5 Violence Against Civilians

set.seed(123456)
violence_data_sf <- grped_acled_sf %>% filter(event_type == "Violence against civilians")

# Initialize an empty list to store the results
violence_G_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(violence_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- violence_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[mmr_owin]
  G_CK.csr <- envelope(filtered_ppp, Gest, nsim = 99)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  violence_G_CK_list[[name]] <- G_CK.csr
}
write_rds(violence_G_CK_list, "data/rds/violence_G_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(violence_G_CK_list)) {
  g_ck <- violence_G_CK_list[[i]]
  plot(g_ck, xlim=c(0,10000), main=names(violence_G_CK_list)[i])
}

5.2.6 Insights

• Majority of the observed G-function plotted for each quarter for the different events lies above the envelope, suggesting that the points, in this case the conflicts, are more clustered than expected under CSR. This means that the distribution of the conflicts is not at random and tends to cluster in specific regions.

• The slope of the observed G-function also generally shows a smooth and constant trend.

• The only exception to the observed trends is the plots for 2021 Q1, which display a more irregular pattern. This irregularity can be attributed to having fewer data points during this time, as the civil war had only just begun in February 2021. When the number of events is smaller, the statistical power of the analysis decreases, leading to less smooth or more erratic G-function patterns. As a result, the G-function for this early quarter may not accurately reflect clustering or dispersion. Given the sparse data from this period, we will not be placing significant weight on the 2021 Q1 results when making decisions about accepting or rejecting the hypothesis.

Conclusion:

• Based on the analysis, the majority of the observed G-function lies outside the envelope, indicating that the conflict events in Myanmar exhibit a pattern of clustering rather than random distribution. This leads us to reject the null hypothesis that the conflicts are randomly distributed. Instead, the results suggest that conflict events are spatially concentrated in specific regions.

5.3 Hotspot of Conflict

5.3.1 Identifying the Hotspot of Conflict

While the KDE map shows a high concentration of conflict in the region where Sagaing, Magway, and Mandalay intersect, it’s challenging to visually pinpoint the exact region with the highest density of conflict events. To overcome this, we employ a quantitative approach using the following code chunk to calculate the number of conflict events within each region.

By counting the number of data points within each region, we observe that the majority of conflict events occur in the Sagaing region. This result aligns with the spatial trends observed in the KDE map, confirming Sagaing as the area with the highest concentration of conflict.

# Create an empty dataframe to store results
conflict_counts <- data.frame(Region = character(), Conflict_Count = integer(), stringsAsFactors = FALSE)

# Get a list of unique regions
regions <- unique(mmrsr_sf$ST)

# Loop through each region
for (region in regions) {
  # Filter the region geometry
  region_sf <- mmrsr_sf %>% filter(ST == region)
  
  # Perform spatial intersection with conflict points
  no_of_conflict_within <- st_intersection(acled_sf, region_sf)
  
  # Count the number of conflicts in this region
  conflict_count <- nrow(no_of_conflict_within)
  
  # Store the results in the dataframe
  conflict_counts <- rbind(conflict_counts, data.frame(Region = region, Conflict_Count = conflict_count))
}

# Arrange the results in descending order of conflict count
conflict_counts_ordered <- conflict_counts %>%
  arrange(desc(Conflict_Count))

write_rds(conflict_counts_ordered, "data/rds/conflict_counts_ordered.rds")

conflict_counts_ordered <- read_rds( "data/rds/conflict_counts_ordered.rds")
print(conflict_counts_ordered)

         Region Conflict_Count
1       Sagaing          11128
2        Magway           4179
3      Mandalay           3603
4  Shan (North)           2938
5        Kachin           2776
6        Yangon           2608
7       Rakhine           2277
8   Tanintharyi           2240
9         Kayin           1817
10          Mon           1677
11         Chin           1557
12 Shan (South)           1430
13        Kayah           1327
14  Bago (East)           1229
15   Ayeyarwady            852
16  Bago (West)            639
17  Nay Pyi Taw            268
18  Shan (East)             56

5.3.2 Sagaing Region

Given that conflict events in Sagaing account for around 25% of the total, we aim to conduct a more detailed analysis to understand the spatial patterns within this significant region. The first step involves isolating the conflict events specific to Sagaing by filtering the data to focus exclusively on this area. For this, instead of the G-function previously used, we will be using the K-function instead since we are now working with a smaller dataset.

5.3.2.1 Overview of Sagaing Region

Since we are focusing solely on Sagaing, here is a more granular map from MIMU which divides Sagaing into districts so we get a better overview of it.

sagaing_map_sf <- mmrsr_district_sf %>% filter(ST == "Sagaing")
tm_shape(sagaing_map_sf) +
  tm_polygons(col = "DT", palette = colors) +  # Apply color palette to polygons
  tm_text("DT", size = 1, col = "black", bg.color = "white", just = c("center", "center"),  xmod = 0, ymod = 0) + tm_layout(main.title = "Sagaing",
            main.title.position = "center",
            main.title.size = 1.6,
            legend.outside = TRUE,
            frame = TRUE)+
    tm_legend(title = "Districts")  # Set custom legend title

This chunk of code creates a Sagaing owin object, converts acled_sf into a ppp object and plots a map displaying only conflict events in Sagaing Region.

sagaing_owin <- as.owin(sagaing_map_sf)
acled_ppp <- as.ppp(acled_sf)
sagaing_ppp <- acled_ppp[sagaing_owin]
plot(sagaing_ppp, pch = 16, cex = 0.5)

5.3.2.2 Hypothesis

Similar to the one we did for the whole of myanmar, we will conduct a hypothesis test using the K-function along with Monte Carlo simulation envelopes. The test will be designed as follows:

• H₀: The distribution of conflict events in Sagaing is randomly distributed (follows Complete Spatial Randomness, CSR).

• H₁: The distribution of conflict events in Sagaing is not randomly distributed.

• If the observed G-function lies outside the envelope, we will reject the null hypothesis, suggesting the presence of spatial clustering or dispersion.

• A 95% confidence interval will be used for the envelopes.

5.3.2.3 Battles

set.seed(123456)
battle_data_sf <- grped_acled_sf %>% filter(event_type == "Battles")

# Initialize an empty list to store the results
sagaing_battle_K_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(battle_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- battle_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[sagaing_owin]
  K_ck.csr <- envelope(filtered_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  sagaing_battle_K_CK_list[[name]] <- K_ck.csr
}

write_rds(sagaing_battle_K_CK_list, "data/rds/sagaing_battle_K_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(sagaing_battle_K_CK_list)) {
  k_ck <- sagaing_battle_K_CK_list[[i]]
  plot(k_ck, . - r ~ r, xlab="d", ylab="k(d)-r", main=names(sagaing_battle_K_CK_list)[i])
}

5.3.2.4 Explosions/ Remote Violence

set.seed(123456)
explosions_data_sf <- grped_acled_sf %>% filter(event_type == "Explosions/Remote violence")

# Initialize an empty list to store the results
sagaing_explosions_K_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(explosions_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- explosions_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[sagaing_owin]
  K_ck.csr <- envelope(filtered_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  sagaing_explosions_K_CK_list[[name]] <- K_ck.csr
}

write_rds(sagaing_explosions_K_CK_list, "data/rds/sagaing_explosions_K_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(sagaing_explosions_K_CK_list)) {
  k_ck <- sagaing_explosions_K_CK_list[[i]]
  plot(k_ck, . - r ~ r, xlab="d", ylab="k(d)-r", main=names(sagaing_explosions_K_CK_list)[i])
}

5.3.2.5 Strategic Developments

set.seed(123456)
strategic_data_sf <- grped_acled_sf %>% filter(event_type == "Strategic developments")

# Initialize an empty list to store the results
sagaing_strategic_K_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(strategic_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- strategic_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[sagaing_owin]
  K_ck.csr <- envelope(filtered_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  sagaing_strategic_K_CK_list[[name]] <- K_ck.csr
}

write_rds(sagaing_strategic_K_CK_list, "data/rds/sagaing_strategic_K_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(sagaing_strategic_K_CK_list)) {
  k_ck <- sagaing_strategic_K_CK_list[[i]]
  plot(k_ck, . - r ~ r, xlab="d", ylab="k(d)-r", main=names(sagaing_strategic_K_CK_list)[i])
}

5.3.2.6 Violence Against Civilians

set.seed(123456)
violence_data_sf <- grped_acled_sf %>% filter(event_type == "Violence against civilians")

# Initialize an empty list to store the results
sagaing_violence_K_CK_list <- list()

# Loop through each combination of year and quarter
for (i in seq_len(nrow(violence_data_sf))) {
  
  # Filter the data for the current iteration
  filtered_sf <- violence_data_sf[i, ]
  
  # Extract the current year, quarter, and event type
  current_year <- filtered_sf$year
  current_quarter <- filtered_sf$quarter
  
  # Convert to ppp object
  filtered_ppp <- as.ppp(st_coordinates(filtered_sf), st_bbox(filtered_sf))
  
  if(any(duplicated(filtered_ppp))){
    filtered_ppp <- rjitter(filtered_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  filtered_ppp <- filtered_ppp[sagaing_owin]
  K_ck.csr <- envelope(filtered_ppp, Kest, nsim = 99, rank = 1, glocal=TRUE)
  
  # Save the result in the list with a unique name
  name <- paste0(current_year, "_Q", current_quarter)
  sagaing_violence_K_CK_list[[name]] <- K_ck.csr
}

write_rds(sagaing_violence_K_CK_list, "data/rds/sagaing_violence_K_CK_list.rds")

par(mar = c(2,1,2,1))
par(mfrow=c(4,4))
# Plotting
for(i in seq_along(sagaing_violence_K_CK_list)) {
  k_ck <- sagaing_violence_K_CK_list[[i]]
  plot(k_ck, . - r ~ r, xlab="d", ylab="k(d)-r", main=names(sagaing_violence_K_CK_list)[i])
}

5.3.2.7 Insights

• Similar to the G-function plots for Myanmar, we can see that most of the plots are identical with majority of the observed K-function being outside and above the envelope with constant and smooth trends.

• The smooth trend of the K-function across the plots implies that the clustering of events is not only prevalent but also persistent over time.

• The only exception to this is yet again 2021 q1, this time we can also observe a stepped pattern with the plots which again is likely attributed to the small amount of data points during this period of time, and hence we will not be placing significant weight on the 2021 Q1 results when making decisions about accepting or rejecting the hypothesis.

• Since the observed K-function are all above the envelope, we can reject the null hypothesis that conflict events in Sagaing is randomly distributed and instead conclude that they are clustered at specific regions. Observing the point pattern map we previously plotted, we can see that most conflict events occur in the following districts: Kanbalu, Shwebo, Sagaing (District, not region), Monywa, Yinmarbi.

6. Quarterly Spatio-Temporal KDE Layers

In this section, we will be conducting spatiotemporal Kernel Density Estimation (KDE). This approach allows us to visualize and analyze the density of conflicts over both space and time, providing a comprehensive view of how conflict activity evolves. By applying KDE to the dataset spanning from 2021 to 2024, we can identify hotspots and trends in conflict intensity across different quarters and regions.

6.1 Battles

6.1.1 Overview

This chunk of code utilizes spattemp.density() to compute the spatiotemporal KDE for all 14 quarters from 2021 to 2024 Q2.

  par(mar = c(2,0,2,0))
  par(mfrow=c(4,4))
  current_event = "Battles"
  
  # Filter the data
  battle_data_sf <- grped_acled_sf %>% filter(event_type == current_event)
  
  # Define time slices for plotting
  tims <- c(
  "2021 Q1", "2021 Q2", "2021 Q3", "2021 Q4",
  "2022 Q1", "2022 Q2", "2022 Q3", "2022 Q4",
  "2023 Q1", "2023 Q2", "2023 Q3", "2023 Q4",
  "2024 Q1", "2024 Q2"
  )
  
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(battle_data_sf), st_bbox(battle_data_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]

  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(tims[x]))
}

6.1.2 By Year

  current_event = "Battles"
  quarter_2021_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2021)
  quarter_2022_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2022)
  quarter_2023_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2023)
  quarter_2024_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2024)

2021

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2021
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2021_sf), st_bbox(quarter_2021_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2022

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2022
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2022_sf), st_bbox(quarter_2022_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2023

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2023
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2023_sf), st_bbox(quarter_2023_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2024

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2024
  
  tims <- c("Q1", "Q2")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2024_sf), st_bbox(quarter_2024_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

6.2 Explosions/Remote violence

6.2.1 Overview

  par(mar = c(2,0,2,0))
  par(mfrow=c(4,4))
  current_event = "Explosions/Remote violence"

  # Filter the data
  explosions_data_sf <- grped_acled_sf %>% filter(event_type == current_event)
  
  # Define time slices for plotting
  tims <- c(
  "2021 Q1", "2021 Q2", "2021 Q3", "2021 Q4",
  "2022 Q1", "2022 Q2", "2022 Q3", "2022 Q4",
  "2023 Q1", "2023 Q2", "2023 Q3", "2023 Q4",
  "2024 Q1", "2024 Q2"
)
  
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(explosions_data_sf), st_bbox(explosions_data_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(tims[x]))
}

6.2.2 By Year

  current_event = "Explosions/Remote violence"
  quarter_2021_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2021)
  quarter_2022_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2022)
  quarter_2023_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2023)
  quarter_2024_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2024)

2021

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2021
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2021_sf), st_bbox(quarter_2021_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2022

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2022
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2022_sf), st_bbox(quarter_2022_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2023

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2023
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2023_sf), st_bbox(quarter_2023_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2024

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2024
  
  tims <- c("Q1", "Q2")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2024_sf), st_bbox(quarter_2024_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

6.3 Strategic developments

6.3.1 Overview

  par(mar = c(2,0,2,0))
  par(mfrow=c(4,4))
  current_event = "Strategic developments"
  
  # Filter the data
  strategic_data_sf <- grped_acled_sf %>% filter(event_type == current_event)
  
  # Define time slices for plotting
  tims <- c(
  "2021 Q1", "2021 Q2", "2021 Q3", "2021 Q4",
  "2022 Q1", "2022 Q2", "2022 Q3", "2022 Q4",
  "2023 Q1", "2023 Q2", "2023 Q3", "2023 Q4",
  "2024 Q1", "2024 Q2"
)
  
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(strategic_data_sf), st_bbox(strategic_data_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(tims[x]))
}

6.3.2 By Year

  current_event = "Strategic developments"
  quarter_2021_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2021)
  quarter_2022_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2022)
  quarter_2023_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2023)
  quarter_2024_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2024)

2021

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2021
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2021_sf), st_bbox(quarter_2021_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2022

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2022
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2022_sf), st_bbox(quarter_2022_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2023

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2023
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2023_sf), st_bbox(quarter_2023_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2024

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2024
  
  tims <- c("Q1", "Q2")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2024_sf), st_bbox(quarter_2024_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

6.4 Violence against civilians

6.4.1 Overview

  par(mar = c(2,0,2,0))
  par(mfrow=c(4,4))
  current_event = "Violence against civilians"
  # Filter the data
  violence_data_sf <- grped_acled_sf %>% filter(event_type == current_event)
  
  # Define time slices for plotting
  tims <- c(
  "2021 Q1", "2021 Q2", "2021 Q3", "2021 Q4",
  "2022 Q1", "2022 Q2", "2022 Q3", "2022 Q4",
  "2023 Q1", "2023 Q2", "2023 Q3", "2023 Q4",
  "2024 Q1", "2024 Q2"
)
  
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(violence_data_sf), st_bbox(violence_data_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(tims[x]))
}

6.4.2 By Year

  current_event = "Violence against civilians"
  quarter_2021_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2021)
  quarter_2022_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2022)
  quarter_2023_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2023)
  quarter_2024_sf <- grped_acled_sf %>% filter(event_type == current_event, year == 2024)

2021

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2021
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2021_sf), st_bbox(quarter_2021_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2022

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2022
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2022_sf), st_bbox(quarter_2022_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2023

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2023
  
  tims <- c("Q1", "Q2", "Q3", "Q4")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2023_sf), st_bbox(quarter_2023_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

2024

  par(mar = c(2,0,2,0))
  par(mfrow=c(2,2))
  year = 2024
  
  tims <- c("Q1", "Q2")
  # Convert to ppp object
  spattemp_ppp <- as.ppp(st_coordinates(quarter_2024_sf), st_bbox(quarter_2024_sf))
  # Remove Duplicates
  if(any(duplicated(spattemp_ppp))){
    rjitter(spattemp_ppp, retry=TRUE, nsim=1, drop=TRUE)
  }
  
  # Convert to owin
  spattemp_ppp <- spattemp_ppp[mmr_owin]
  
  # Compute the spatiotemporal Kernel Density Estimate (KDE)
  st_kde <- spattemp.density(spattemp_ppp)
  for(x in seq_along(tims)){
    plot(st_kde, x, 
     override.par=FALSE, 
     fix.range=TRUE, 
     main=paste(year,tims[x]))}

7. 2nd-Order Spatio-temporal Point Patterns Analysis

The specific questions we aim to address are whether the locations of conflict in Myanmar are spatio-temporally independent. Using second-order spatio-temporal point pattern analysis, we will examine if the conflicts exhibit significant clustering or dispersion across geographic areas and over time.

This section is definitely the toughest for me as it was hard to find resources and documentations online to perform the relevant analysis.

7.1 Learning Process

7.1.1 Starting Out

To start off, I researched on the potential Second Order Spatio-temporal Point Patterns Analysis methods I could use and ended up shortlisting out these options:

• Spatio-temporal K-function: This function extends the traditional K-function to consider both space and time, allowing us to assess how conflicts cluster over time and across different locations.

• Ripley’s K-function: Specifically, the K-function can be adapted for spatio-temporal analysis to evaluate the density of events in both dimensions. This helps in identifying patterns of clustering or dispersion over specified time intervals.

• Log-Gaussian Cox Process (LGCP): This is a flexible model that accounts for spatial dependence in point patterns and allows for varying intensity across space and time. LGCP can help in modeling the intensity of conflict events as a function of spatial covariates, providing insights into the underlying processes driving the distribution of conflicts.

Instinctively, I was inclined to stick with the K-function since it is a method I am familiar with and have applied before, making it a comfortable starting point for my analysis. However, I also recognized that the term “Log-Gaussian Cox Process” appeared frequently in the literature I reviewed, which prompted me to include it in my shortlist.

7.1.2 Failed Attempts

Despite trying many different functions and methods, not many of them came through. A method that I had initially thought worked with the help of ChatGPT was using time (in this case quarter number) to mark the point patterns. However, I soon realised after that the standard Kest function that we used in class does not seem to take into account marks.

Using Mark for K-Function

current_event <- "Battles"

quarter_to_continuous <- function(year, quarter) {
  # Calculate the number of quarters since a reference point (e.g., 2000 Q1)
  reference_year <- 2021
  quarters_since_reference <- (year - reference_year) * 4 + quarter
  return(quarters_since_reference)
}

grped_acled_quartered_sf <- acled_sf %>%   mutate(
                                                          year = as.numeric(year),
                                                          quarter = as.numeric(quarter),
                                                          t = quarter_to_continuous(year, quarter)
                                                        )
grped_acled_quartered_sf <- grped_acled_quartered_sf %>% filter(event_type == current_event)

# Step 2: Convert the data to a spatio-temporal point pattern
# Define the spatial window (using your bounding box)

# Create a spatio-temporal point pattern (ppp is the spatial part, time is added separately)
ppp_data <- as.ppp(st_coordinates(grped_acled_quartered_sf), W = mmr_owin)
marks(ppp_data) <- data.frame(time = acled_sf$t)

Kst_result <- Kest(ppp_data, tmax = 1, rmax = 10000)

# Step 4: Visualize and interpret the results
plot(Kst_result)

7.2 Space-Time Inhomogeneous K-function

The method that I eventually ended up going for, and succeeding in getting to actually work (or at least have some semblance of success) is using STIKhat() from the package stpp which is said to be designed for “Space-Time Point Pattern Simulation, Visualisation and Analysis”.

STIKhat(): Compute an estimate of the Space-Time Inhomogeneous K-function.

This chunk of code converts acled data to continuous quarter for time slices and gets rid of duplicate rows with the same event type, geometry and t (overall quarter no, e.g. 2024 Q2 –> 14).

  quarter_to_continuous <- function(year, quarter) {
    # Calculate the number of quarters since a reference point (e.g., 2000 Q1)
    reference_year <- 2021
    quarters_since_reference <- (year - reference_year) * 4 + quarter
    return(quarters_since_reference)
  }
  
  grped_acled_quartered_sf <- acled_sf %>%   mutate(year = as.numeric(year),
                                                    quarter = as.numeric(quarter),
                                                    t = quarter_to_continuous(year, quarter)
                                                    ) %>% distinct(event_type, geometry, t, .keep_all = TRUE)

Important

There seems to be an issue where the plots appear in the viewer of Rstudio but yet I am unable to have it rendered and appear on the Quarto Document. Hence I have made the choice to instead download and display it instead.

This first chunk of code utilises STIKhat() from the stpp package to compute an estimate of the Space-Time Inhomogeneous K-function and plotK() to plot the graph.

The second chunk of code then plots the different graphical representation of the estimate of the Space-Time Inhomogeneous K-function: Projection, Perspective, Image so that we are able to get a view of the different graphs.

Battles

  current_event <- "Battles"

  battle_quartered_sf <- grped_acled_quartered_sf %>% filter(event_type == current_event)
  coords_matrix <- st_coordinates(battle_quartered_sf)
  x_coords <- coords_matrix[, "X"]
  y_coords <- coords_matrix[, "Y"]
  
  stpp_data <- as.3dpoints(x_coords,y_coords, battle_quartered_sf$t)
  
  Kst_result <- stpp::STIKhat(xyt=stpp_data)

  write_rds(Kst_result, "data/rds/battle_stik_result.rds")

plotK(battle_stik_result)
plotK(battle_stik_result, type ="persp")
plotK(battle_stik_result, type ="image")

Battle Projection

Battle Perspective

Battle Image

Explosions/ Remote Violence

  current_event <- "Explosions/Remote violence"

  explosions_quartered_sf <- grped_acled_quartered_sf %>% filter(event_type == current_event)
  
  coords_matrix <- st_coordinates(explosions_quartered_sf)
  x_coords <- coords_matrix[, "X"]
  y_coords <- coords_matrix[, "Y"]
  
  stpp_data <- as.3dpoints(x_coords,y_coords, explosions_quartered_sf$t)
  
  Kst_result <- stpp::STIKhat(xyt=stpp_data)

  write_rds(Kst_result, "data/rds/explosions_stik_result.rds")

plotK(explosions_stik_result)
plotK(explosions_stik_result, type ="persp")
plotK(explosions_stik_result, type ="image")

Explosions Projection

Explosions Perspective

Explosions Image

Strategic Developments

  current_event <- "Strategic developments"

  strategic_quartered_sf <- grped_acled_quartered_sf %>% filter(event_type == current_event)
  
  coords_matrix <- st_coordinates(strategic_quartered_sf)
  x_coords <- coords_matrix[, "X"]
  y_coords <- coords_matrix[, "Y"]
  
  stpp_data <- as.3dpoints(x_coords,y_coords, strategic_quartered_sf$t)
  
  Kst_result <- stpp::STIKhat(xyt=stpp_data)

  write_rds(Kst_result, "data/rds/strategic_stik_result.rds")

plotK(strategic_stik_result)
plotK(strategic_stik_result, type ="persp")
plotK(strategic_stik_result, type ="image")

Strategic Projection

Strategic Perspective

Strategic Image

Violence Against Civilians

  current_event <- "Violence against civilians"
  violence_quartered_sf <- grped_acled_quartered_sf %>% filter(event_type == current_event)
  
  coords_matrix <- st_coordinates(violence_quartered_sf)
  x_coords <- coords_matrix[, "X"]
  y_coords <- coords_matrix[, "Y"]
  
  stpp_data <- as.3dpoints(x_coords,y_coords, violence_quartered_sf$t)
  
  Kst_result <- stpp::STIKhat(xyt=stpp_data)

  write_rds(Kst_result, "data/rds/violence_stik_result.rds")

plotK(violence_stik_result)
plotK(violence_stik_result, type ="persp")
plotK(violence_stik_result, type ="image")

Violence Projection

Violence Perspective

Violence Image

7.3 Interpreting the Plots

8. OpenStreetmap

9. Overall Insights

10. Reflections

The biggest challenge I faced in this assignment was definitely in the interpretation of the plots and the analysis. It was difficult to find online resources on how to interpret it, and while there were numerous technical guides on how to calculate spatial statistics, very few focused on providing clear explanations for what the results actually mean in a real-world context.

Understanding how to interpret functions like the G-function and the significance of clustering versus dispersion required a deep dive into spatial theory, which was not always straightforward.

Another layer of complexity was in distinguishing between statistical significance and practical relevance. While my results showed clear clustering, translating that into meaningful insights—such as understanding why certain areas experience more conflicts—was not immediately obvious. This challenge highlighted the importance of combining spatial analysis with domain knowledge, in this case, political or socio-economic factors driving conflict.

Furthermore, I also had to carefully consider the assumptions behind the models I used, particularly the assumption of Complete Spatial Randomness (CSR). In reality, conflicts are rarely random, and this led to questions about how appropriate these models are for analyzing such phenomena. It pushed me to think critically about the limitations of the tools and methods and reminded me that spatial analysis is as much about the context and the data as it is about the statistical techniques.

In the end, this assignment taught me the importance of not just relying on statistical results but also developing the ability to interpret and contextualize those results in light of the real-world phenomena they are meant to represent. It also reinforced the need for interdisciplinary approaches, combining quantitative analysis with qualitative insights.

11. References

Baddeley, A., Rubak, E., & Turner, R. (2016). Spatial point patterns: Methodology and applications with R. CHAPMAN & HALL CRC.